
rm(list=ls())

  
  library(ggplot2)
  library(gridExtra)
  library(haven)
  library(margins)
  library(tidyverse)
  library(ggridges)
  library(magrittr)
  library(broom)
  library(rddtools)
  library(rdrobust)
  library(devtools)
  library(patchwork)
  library(stargazer)
  library(VGAM)
  library(ggeffects)
library(stats)
library(Cairo)
  
###################################################################Turnout 2012############################################################################


##### Inferential Analyses #####
setwd("E:\\Dropbox\\RR Discrim Turnout\\REPLICATION FOLDER") # REPLACE WITH DIRECTORY TO PROJECT

ANES2012 <- read_dta("ANES 2012Rec.dta")
whites12 <- subset(ANES2012, white==1)

interactiondemcontrols12 <- glm(voted ~ reldep * rr + age + income + education + female + resimobi + ownhome, data = whites12, family = "binomial")
summary(interactiondemcontrols12)

interactionfullcontrols12 <- glm(voted ~ reldep * rr + age + income + education + female + resimobi + ownhome + strpty + church_attoft + union + contacted , data = whites12, family = "binomial")
summary(interactionfullcontrols12)

test <- with(whites12, data.frame(rr=c(0,1,0,1),
                                  reldep=c(0,0,1,1),
                                  age=mean(age,na.rm=TRUE),
                                  income=mean(income,na.rm=TRUE),
                                  education=mean(education,na.rm=TRUE),
                                  female=mean(female,na.rm=TRUE),
                                  resimobi=mean(resimobi,na.rm=TRUE),
                                  ownhome=mean(ownhome,na.rm=TRUE),
                                  strpty=mean(strpty,na.rm=TRUE),
                                 church_attoft=mean(church_attoft,na.rm=TRUE),
                                  union=mean(union,na.rm=TRUE),
                                  contacted=mean(contacted,na.rm=TRUE)))

predict(interactionfullcontrols12, test, type="response")

###################################################################Turnout 2016############################################################################



ANES2016 <- read_dta("ANES2016Rec.dta")
whites16 <- subset(ANES2016, white==1)

interactiondemcontrols16 <- glm(voted ~ reldep * rr + age + income + education + female + resimobi + ownhome, data = whites16, family = "binomial")
summary(interactiondemcontrols16)

interactionfullcontrols16 <- glm(voted ~ reldep * rr + age + income + education + female + resimobi + ownhome + strpty + church_attoft + union + contacted , data = whites16, family = "binomial")
summary(interactionfullcontrols16)

test <- with(whites16, data.frame(rr=c(0,1,0,1),
                                  reldep=c(0,0,1,1),
                                  age=mean(age,na.rm=TRUE),
                                  income=mean(income,na.rm=TRUE),
                                  education=mean(education,na.rm=TRUE),
                                  female=mean(female,na.rm=TRUE),
                                  resimobi=mean(resimobi,na.rm=TRUE),
                                  ownhome=mean(ownhome,na.rm=TRUE),
                                  strpty=mean(strpty,na.rm=TRUE),
                                  church_attoft=mean(church_attoft,na.rm=TRUE),
                                  union=mean(union,na.rm=TRUE),
                                  contacted=mean(contacted,na.rm=TRUE)))

predict(interactionfullcontrols16, test, type="response")

###################################################################Turnout 2020############################################################################

ANES2020 <- read_dta("ANES 2020 Rec.dta")
whites20 <- subset(ANES2020, r_white==1)

interactiondemcontrols20 <- glm(voted ~ reldep * rr + age + income + education + female + resimobi + ownhome, data = whites20, family = "binomial")
summary(interactiondemcontrols20)

interactionfullcontrols20 <- glm(voted ~ reldep * rr + age + income + education + female + resimobi + ownhome + strpty + church_attoft + union + contacted , data = whites20, family = "binomial")
summary(interactionfullcontrols20)

test <- with(whites20, data.frame(rr=c(0,1,0,1),
                                  reldep=c(0,0,1,1),
                                  age=mean(age,na.rm=TRUE),
                                  income=mean(income,na.rm=TRUE),
                                  education=mean(education,na.rm=TRUE),
                                  female=mean(female,na.rm=TRUE),
                                  resimobi=mean(resimobi,na.rm=TRUE),
                                  ownhome=mean(ownhome,na.rm=TRUE),
                                  strpty=mean(strpty,na.rm=TRUE),
                                  church_attoft=mean(church_attoft,na.rm=TRUE),
                                  union=mean(union,na.rm=TRUE),
                                  contacted=mean(contacted,na.rm=TRUE)))

predict(interactionfullcontrols20, test, type="response")

###################################################################Participation 2012############################################################################

eb_interactiondemcontrols12 <- lm(eb_index7 ~ reldep * rr + age + income + education + female + resimobi + ownhome, data = whites12)
summary(eb_interactiondemcontrols12)

eb_interactionfullcontrols12 <- lm(eb_index7 ~ reldep * rr + age + income + education + female + resimobi + ownhome + strpty + church_attoft + union + contacted, data = whites12)
summary(eb_interactionfullcontrols12)

test <- with(whites12, data.frame(rr=c(0,1,0,1),
                                  reldep=c(0,0,1,1),
                                  age=mean(age,na.rm=TRUE),
                                  income=mean(income,na.rm=TRUE),
                                  education=mean(education,na.rm=TRUE),
                                  female=mean(female,na.rm=TRUE),
                                  resimobi=mean(resimobi,na.rm=TRUE),
                                  ownhome=mean(ownhome,na.rm=TRUE),
                                  strpty=mean(strpty,na.rm=TRUE),
                                  church_attoft=mean(church_attoft,na.rm=TRUE),
                                  union=mean(union,na.rm=TRUE),
                                  contacted=mean(contacted,na.rm=TRUE)))

predict(eb_interactionfullcontrols12, test, type="response")*7

ce_interactiondemcontrols12 <- lm(ce_index ~ reldep * rr + age + income + education + female + resimobi + ownhome, data = whites12)
summary(ce_interactiondemcontrols12)

ce_interactionfullcontrols12 <- lm(ce_index ~ reldep * rr + age + income + education + female + resimobi + ownhome + strpty + church_attoft + union + contacted, data = whites12)
summary(ce_interactionfullcontrols12)

test <- with(whites12, data.frame(rr=c(0,1,0,1),
                                  reldep=c(0,0,1,1),
                                  age=mean(age,na.rm=TRUE),
                                  income=mean(income,na.rm=TRUE),
                                  education=mean(education,na.rm=TRUE),
                                  female=mean(female,na.rm=TRUE),
                                  resimobi=mean(resimobi,na.rm=TRUE),
                                  ownhome=mean(ownhome,na.rm=TRUE),
                                  strpty=mean(strpty,na.rm=TRUE),
                                  church_attoft=mean(church_attoft,na.rm=TRUE),
                                  union=mean(union,na.rm=TRUE),
                                  contacted=mean(contacted,na.rm=TRUE)))

predict(ce_interactionfullcontrols12, test, type="response")

###################################################################Participation 2016############################################################################

eb_interactiondemcontrols16 <- lm(eb_index7 ~ reldep * rr + age + income + education + female + resimobi + ownhome, data = whites16)
summary(eb_interactiondemcontrols16)

eb_interactionfullcontrols16 <- lm(eb_index7 ~ reldep * rr + age + income + education + female + resimobi + ownhome + strpty + church_attoft + union + contacted, data = whites16)
summary(eb_interactionfullcontrols16)

test <- with(whites16, data.frame(rr=c(0,1,0,1),
                                  reldep=c(0,0,1,1),
                                  age=mean(age,na.rm=TRUE),
                                  income=mean(income,na.rm=TRUE),
                                  education=mean(education,na.rm=TRUE),
                                  female=mean(female,na.rm=TRUE),
                                  resimobi=mean(resimobi,na.rm=TRUE),
                                  ownhome=mean(ownhome,na.rm=TRUE),
                                  strpty=mean(strpty,na.rm=TRUE),
                                  church_attoft=mean(church_attoft,na.rm=TRUE),
                                  union=mean(union,na.rm=TRUE),
                                  contacted=mean(contacted,na.rm=TRUE)))

predict(eb_interactionfullcontrols16, test, type="response")*8

ce_interactiondemcontrols16 <- lm(ce_index ~ reldep * rr + age + income + education + female + resimobi + ownhome, data = whites16)
summary(ce_interactiondemcontrols16)

ce_interactionfullcontrols16 <- lm(ce_index ~ reldep * rr + age + income + education + female + resimobi + ownhome + strpty + church_attoft + union + contacted, data = whites16)
summary(ce_interactionfullcontrols16)

test <- with(whites16, data.frame(rr=c(0,1,0,1),
                                  reldep=c(0,0,1,1),
                                  age=mean(age,na.rm=TRUE),
                                  income=mean(income,na.rm=TRUE),
                                  education=mean(education,na.rm=TRUE),
                                  female=mean(female,na.rm=TRUE),
                                  resimobi=mean(resimobi,na.rm=TRUE),
                                  ownhome=mean(ownhome,na.rm=TRUE),
                                  strpty=mean(strpty,na.rm=TRUE),
                                  church_attoft=mean(church_attoft,na.rm=TRUE),
                                  union=mean(union,na.rm=TRUE),
                                  contacted=mean(contacted,na.rm=TRUE)))

predict(ce_interactionfullcontrols16, test, type="response")*7

###################################################################Participation 2020############################################################################

eb_interactiondemcontrols20 <- lm(eb_index7 ~ reldep * rr + age + income + education + female + resimobi + ownhome, data = whites20)
summary(eb_interactiondemcontrols20)

eb_interactionfullcontrols20 <- lm(eb_index7 ~ reldep * rr + age + income + education + female + resimobi + ownhome + strpty + church_attoft + union + contacted, data = whites20)
summary(eb_interactionfullcontrols20)

test <- with(whites20, data.frame(rr=c(0,1,0,1),
                                  reldep=c(0,0,1,1),
                                  age=mean(age,na.rm=TRUE),
                                  income=mean(income,na.rm=TRUE),
                                  education=mean(education,na.rm=TRUE),
                                  female=mean(female,na.rm=TRUE),
                                  resimobi=mean(resimobi,na.rm=TRUE),
                                  ownhome=mean(ownhome,na.rm=TRUE),
                                  strpty=mean(strpty,na.rm=TRUE),
                                  church_attoft=mean(church_attoft,na.rm=TRUE),
                                  union=mean(union,na.rm=TRUE),
                                  contacted=mean(contacted,na.rm=TRUE)))

predict(eb_interactionfullcontrols20, test, type="response")*8


ce_interactiondemcontrols20 <- lm(ce_index ~ reldep * rr + age + income + education + female + resimobi + ownhome, data = whites20)
summary(ce_interactiondemcontrols20)

ce_interactionfullcontrols20 <- lm(ce_index ~ reldep * rr + age + income + education + female + resimobi + ownhome + strpty + church_attoft + union + contacted, data = whites20)
summary(ce_interactionfullcontrols20)

test <- with(whites20, data.frame(rr=c(0,1,0,1),
                                  reldep=c(0,0,1,1),
                                  age=mean(age,na.rm=TRUE),
                                  income=mean(income,na.rm=TRUE),
                                  education=mean(education,na.rm=TRUE),
                                  female=mean(female,na.rm=TRUE),
                                  resimobi=mean(resimobi,na.rm=TRUE),
                                  ownhome=mean(ownhome,na.rm=TRUE),
                                  strpty=mean(strpty,na.rm=TRUE),
                                  church_attoft=mean(church_attoft,na.rm=TRUE),
                                  union=mean(union,na.rm=TRUE),
                                  contacted=mean(contacted,na.rm=TRUE)))

predict(ce_interactionfullcontrols20, test, type="response")*13


##############################################################Table One###################################################
descriptivesdf <- as.data.frame(cbind(whites12$rr, whites12$reldep, whites16$rr, whites16$reldep, whites20$rr, whites20$reldep))

stargazer(descriptivesdf,
          title="Table 1: Descriptives of Independent Variables",
          align=TRUE,
          covariate.labels = c("Racial Resentment 2012",
                               "Relative Deprivation 2012",
                               "Racial Resentment 2016",
                               "Relative Deprivation 2016",
                               "Racial Resentment 2020",
                               "Relative Deprivation 2020"),
          out = "Table B1.html")


descriptivesdf <- as.data.frame(cbind(whites12$voted,
                                      whites12$eb_index7,
                                      whites12$ce_index,
                                      whites16$voted,
                                      whites16$eb_index7,
                                      whites16$ce_index,
                                      whites20$voted,
                                      whites20$eb_index7,
                                      whites20$ce_index))

stargazer(descriptivesdf,
          title="Table 2: Descriptives of Independent Variables",
        align=TRUE,
          covariate.labels = c("Voted 2012",
                               "Political Participation 2012",
                               "Civic Engagement 2012",
                               "Voted 2016",
                               "Political Participation 2016",
                              "Civic Engagement 2016",
                               "Voted 2020",
                               "Political Participation 2020",
                              "Civic Engagement 2020"),
         out = "Table B2.html")

library(psych)

psych::alpha(whites12[,c("eb_contrep",
                         "eb_march",
                         "eb_contrmoney",
                         "eb_campaign",
                         "eb_rally",
                         "eb_talk",
                         "eb_sign")]) # .66

psych::alpha(whites12[,c("ce_letter",
                         "ce_message",
                         "ce_call",
                         "ce_controrg",
                         "ce_contrrelig",
                         "ce_board",
                         "ce_petition1",
                         "ce_petition2",
                         "eb_march",
                         "eb_contrep")]) # .67

psych::alpha(na.omit(whites16[,c("eb_talk",
                         "eb_protest",
                         "eb_campaignlit",
                         "eb_moneycand",
                         "eb_moneyparty",
                         "eb_contactr",
                         "registeredsum",
                         "voted")])) # .58

psych::alpha(whites16[,c("v162195_R",
                         "v162196_R",
                         "v162197_R",
                         "v162198_R",
                         "v162200_R",
                         "v162202_R",
                         "v162204_R")]) #.72

psych::alpha(whites20[,c("V202009_R",
                         "V202013_R",
                         "V202014_R",
                         "V202015_R",
                         "V202016_R",
                         "V202017_R",
                         "V202019_R",
                         "V202021_R")])



psych::alpha(whites20[,c("V202024_R",
                         "V202025_R",
                         "V202026_R",
                         "V202027_R",
                         "V202028_R",
                         "V202029_R",
                         "V202030_R",
                         "V202031_R",
                         "V202032_R",
                         "V202033_R",
                         "V202034_R",
                         "V202036_R",
                         "V202023_R",
                         "V202040_R")]) #.73



##############################################################FIGURE ONE############################################################################

cor(whites12$rr, whites12$reldep, use = "pairwise.complete.obs")

TWELVE <- (ggplot(whites12, aes(rr, reldep)) +
            stat_smooth(size = 1.5, method = "loess", color = "black") +
            xlab("Racial Resentment") +
            ylab("Relative Discrimination") + 
            ylim(0, 1) +             ggtitle("2012") + 
             geom_jitter(alpha=.1) +
             theme_classic() +
            theme(
              plot.title=element_text(hjust=.5, size=14, face="bold")))
  
TWELVE

cor(whites16$rr, whites16$reldep, use = "pairwise.complete.obs")

SIXTEEN <- (ggplot(whites16, aes(rr, reldep)) +
              stat_smooth(size = 1.5, method = "loess", color = "black") +
              xlab("Racial Resentment") +
              ylab("Relative Discrimination") + 
              ylim(0, 1) +             ggtitle("2016") + 
              geom_jitter(alpha=.1) +
              theme_classic() +
              theme(
                plot.title=element_text(hjust=.5, size=14, face="bold")))
SIXTEEN

TWENTY <- (ggplot(whites20, aes(rr, reldep)) +
              stat_smooth(size = 1.5, method = "loess", color = "black") +
              xlab("Racial Resentment") +
              ylab("Relative Discrimination") + 
              ylim(0, 1) +             ggtitle("2020") + 
              geom_jitter(alpha=.05) +
              theme_classic() +
              theme(
                plot.title=element_text(hjust=.5, size=14, face="bold")))
TWENTY

mytable <- 1:6
dim(mytable) <- c(3,2)

mytable[1,1] <- "2012"
mytable[2,1] <- "2016"
mytable[3,1] <- "2020"


mytable[1,2] <- round(cor(whites12$rr, whites12$reldep, use = "pairwise.complete.obs"),3)
mytable[2,2] <- round(cor(whites16$rr, whites16$reldep, use = "pairwise.complete.obs"),3)
mytable[3,2] <- round(cor(whites20$rr, whites20$reldep, use = "pairwise.complete.obs"),3)

colnames(mytable) <- c("Year", "r")

F1 <- (TWELVE + SIXTEEN + TWENTY + gridExtra::tableGrob(mytable)) +
       plot_layout(ncol =2)

png("Figure B1.png", width = 500, height = 500, units = "px")
wrap_elements(F1) 
dev.off()


####################################################Table 1#####################################################################

stargazer(interactiondemcontrols12,
          interactionfullcontrols12,
          interactiondemcontrols16,
          interactionfullcontrols16,
          interactiondemcontrols20,
          interactionfullcontrols20,
          title="Table 1: Racial Resentment and Relative Discrimination Interactively Predict Turnout",
          align=TRUE,
          dep.var.labels.include=FALSE,
          dep.var.caption = "",
          omit.stat = c("ser"),
          order = c("reldep:rr", "rr", "reldep", "age", "income", "education", "female", "resimobi", "ownhome", "strpty", "church_attoft", "union", "contacted", "Constant"),
          covariate.labels = c("Interaction Term",
                               "Rel. Discrimination",
                               "Racial Resentment",
                               "Age",
                               "Income",
                               "Education",
                               "Female",
                               "Resid. Mobil.",
                               "Home Owner",
                               "PID Strength",
                               "Freq. Church Att.",
                               "Union Mem. in HH",
                               "Contacted",
                               "Constant"),
          star.char = c("*"),
          star.cutoffs = c(0.05),
          notes = c("* p<0.05"), 
          notes.append = F,
          out = "Table 1.html")


####################################################Table 2#####################################################################

stargazer(eb_interactiondemcontrols12,
          eb_interactionfullcontrols12,
          eb_interactiondemcontrols16,
          eb_interactionfullcontrols16,
          eb_interactiondemcontrols20,
          eb_interactionfullcontrols20,
          title="Table 2: Racial Resentment and Relative Discrimination Interactively Predict Political Engagement",
          align=TRUE,
          dep.var.labels.include=FALSE,
          dep.var.caption = "",
          omit.stat = c("ser"),
          order = c("reldep:rr", "rr", "reldep", "age", "income", "education", "female", "resimobi", "ownhome", "strpty", "church_attoft", "union", "contacted", "Constant"),
          covariate.labels = c("Interaction Term",
                               "Rel. Discrimination",
                               "Racial Resentment",
                               "Age",
                               "Income",
                               "Education",
                               "Female",
                               "Resid. Mobil.",
                               "Home Owner",
                               "PID Strength",
                               "Freq. Church Att.",
                               "Union Mem. in HH",
                               "Contacted",
                               "Constant"),
          star.char = c("*"),
          star.cutoffs = c(0.05),
          notes = c("* p<0.05"), 
          notes.append = F,
          out = "Table 2.html")


####################################################Table 3#####################################################################

stargazer(ce_interactiondemcontrols12,
          ce_interactionfullcontrols12,
          ce_interactiondemcontrols16,
          ce_interactionfullcontrols16,
          ce_interactiondemcontrols20,
          ce_interactionfullcontrols20,
          title="Table 3: Racial Resentment and Relative Discrimination Interactively Predict Civic Engagement",
          align=TRUE,
          dep.var.labels.include=FALSE,
          dep.var.caption = "",
          omit.stat = c("ser"),
          order = c("reldep:rr", "rr", "reldep", "age", "income", "education", "female", "resimobi", "ownhome", "strpty", "church_attoft", "union", "contacted", "Constant"),
          covariate.labels = c("Interaction Term",
                               "Rel. Discrimination",
                               "Racial Resentment",
                               "Age",
                               "Income",
                               "Education",
                               "Female",
                               "Resid. Mobil.",
                               "Home Owner",
                               "PID Strength",
                               "Freq. Church Att.",
                               "Union Mem. in HH",
                               "Contacted",
                               "Constant"),
          star.char = c("*"),
          star.cutoffs = c(0.05),
          notes = c("* p<0.05"), 
          notes.append = F,
          out = "Table 3.html")




######################################figure 1###################################################

summary(interactionfullcontrols12)
cdat <- cplot(interactionfullcontrols12, "rr", draw = FALSE)
head(cdat)
rr12 <- ggplot(cdat, aes(x = xvals)) + 
  geom_line(aes(y = yvals)) +
  geom_line(aes(y = upper), linetype = 2) +
  geom_line(aes(y = lower), linetype = 2) +
  ylim(0, 1) + 
  ggtitle("Racial Resentment") +
  xlab("Racial Resentment") + ylab("pr(Voted)") +
  theme_classic() +
  theme(
    plot.title=element_text(hjust=.5, size=12, face="bold"),
    axis.text = element_text(color="black"),
    axis.ticks = element_line(color="black"))

summary(interactionfullcontrols12)
cdat <- cplot(interactionfullcontrols12, "reldep", draw = FALSE)
head(cdat)
rd12 <- ggplot(cdat, aes(x = xvals)) + 
  geom_line(aes(y = yvals)) +
  geom_line(aes(y = upper), linetype = 2) +
  geom_line(aes(y = lower), linetype = 2) +
  ylim(0, 1) + 
  ggtitle("Rel. Discrimination") +
  xlab("Rel. Discrimination") + ylab("pr(Voted)") +
  theme_classic() +
  theme(
    plot.title=element_text(hjust=.5, size=12, face="bold"),
    axis.text = element_text(color="black"),
    axis.ticks = element_line(color="black"))




llabs <- c("0" = "Rel. Discrim. = 0",
           "0.25" = "Rel. Discrim. = .25",
           "0.5" = "Rel. Discrim. = .5",
           "0.75" = "Rel. Discrim. = .75",
           "1" = "Rel. Discrim. = 1")

mydf <- ggpredict(interactionfullcontrols12, terms = c("rr", "reldep"), ci.lvl = 0.95)
int12 <- ggplot(mydf, aes(x, predicted, group = group)) +  geom_line() +
  geom_line(aes(y = conf.high), linetype = 2) +
  geom_line(aes(y = conf.low), linetype = 2) + facet_wrap(~group, nrow = 3, labeller=labeller(group = llabs)) +
  xlab("Racial Resentment") + ylab("pr(Voted)") +
  scale_x_continuous(breaks = c(0,
                                .5,
                                1)) +
  ggtitle("Rel. Discrimination *\nRacial Resentment") +
  theme_classic() +
  theme(
    plot.title=element_text(hjust=.5, size=12, face="bold"),
    axis.text = element_text(color="black"),
    axis.ticks = element_line(color="black"))
  
int12
  

MFX12 <- (int12 | (rr12 / rd12)) + plot_layout()
MFX12

png("Figure 1.png", width = 750, height = 500, units = "px")
wrap_elements(MFX12) 
dev.off()

ggsave(wrap_elements(MFX12) ,
       filename = "Figure 1.pdf",
       device=cairo_pdf,
       dpi = 300,
       units = "in",
       width = 6,
       height = 4)

######################################figure 2###################################################

rr16 <- ggpredict(interactionfullcontrols16, terms = c("rr"), ci.lvl = 0.95)
rr16 <- ggplot(rr16, aes(x, predicted)) +  geom_line() +
  geom_line(aes(y = conf.high), linetype = 2) +
  geom_line(aes(y = conf.low), linetype = 2) +
  xlab("Racial Resentment") + ylab("pr(Voted)") +
  scale_x_continuous(breaks = c(0,
                                .5,
                                1)) +
  ggtitle("Racial Resentment") +
  ylim(c(0,1)) +
  theme_classic() +
  theme(
    plot.title=element_text(hjust=.5, size=12, face="bold"),
    axis.text = element_text(color="black"),
    axis.ticks = element_line(color="black"))


rd16 <- ggpredict(interactionfullcontrols16, terms = c("reldep"), ci.lvl = 0.95)
rd16 <- ggplot(rd16, aes(x, predicted)) +  geom_line() +
  geom_line(aes(y = conf.high), linetype = 2) +
  geom_line(aes(y = conf.low), linetype = 2) +
  xlab("Rel. Discrimination") + ylab("pr(Voted)") +
  scale_x_continuous(breaks = c(0,
                                .5,
                                1)) +
  ggtitle("Rel. Discrimination") +
  ylim(c(0,1)) +
  theme_classic() +
  theme(
    plot.title=element_text(hjust=.5, size=12, face="bold"),
    axis.text = element_text(color="black"),
    axis.ticks = element_line(color="black"))


llabs <- c("0" = "Rel. Discrim. = 0",
           "0.333333343267441" = "Rel. Discrim. = .3",
           "0.5" = "Rel. Discrim. = .5",
           "0.666666686534882" = "Rel. Discrim. = .66",
           "1" = "Rel. Discrim. = 1")

mydf <- ggpredict(interactionfullcontrols16, terms = c("rr", "reldep[0, 0.5, 1]"), ci.lvl = 0.95)
int16 <- ggplot(mydf, aes(x, predicted, group = group)) +  geom_line() +
  geom_line(aes(y = conf.high), linetype = 2) +
  geom_line(aes(y = conf.low), linetype = 2) + facet_wrap(~group, nrow = 2, labeller=labeller(group = llabs)) +
  xlab("Racial Resentment") + ylab("pr(Voted)") +
  scale_x_continuous(breaks = c(0,
                                .5,
                                1)) +
  ylim(c(0,1)) +
  ggtitle("Rel. Discrimination *\nRacial Resentment") +
  theme_classic() +
  theme(
    plot.title=element_text(hjust=.5, size=12, face="bold"),
    axis.text = element_text(color="black"),
    axis.ticks = element_line(color="black"))

int16


MFX16 <- (int16 | (rr16 / rd16)) + plot_layout()
MFX16

png("Figure 2.png", width = 750, height = 500, units = "px")
wrap_elements(MFX16) 
dev.off()

ggsave(wrap_elements(MFX16) ,
       filename = "Figure 2.pdf",
       device=cairo_pdf,
       dpi = 300,
       units = "in",
       width = 6,
       height = 4)


######################################figure 3###################################################

rr20 <- ggpredict(interactionfullcontrols20, terms = c("rr"), ci.lvl = 0.95)
rr20 <- ggplot(rr20, aes(x, predicted)) +  geom_line() +
  geom_line(aes(y = conf.high), linetype = 2) +
  geom_line(aes(y = conf.low), linetype = 2) +
  xlab("Racial Resentment") + ylab("pr(Voted)") +
  scale_x_continuous(breaks = c(0,
                                .5,
                                1)) +
  ggtitle("Racial Resentment") +
  ylim(c(0,1)) +
  theme_classic() +
  theme(
    plot.title=element_text(hjust=.5, size=12, face="bold"),
    axis.text = element_text(color="black"),
    axis.ticks = element_line(color="black"))


rd20 <- ggpredict(interactionfullcontrols20, terms = c("reldep"), ci.lvl = 0.95)
rd20 <- ggplot(rd20, aes(x, predicted)) +  geom_line() +
  geom_line(aes(y = conf.high), linetype = 2) +
  geom_line(aes(y = conf.low), linetype = 2) +
  xlab("Rel. Discrimination") + ylab("pr(Voted)") +
  scale_x_continuous(breaks = c(0,
                                .5,
                                1)) +
  ggtitle("Rel. Discrimination") +
  ylim(c(0,1)) +
  theme_classic() +
  theme(
    plot.title=element_text(hjust=.5, size=12, face="bold"),
    axis.text = element_text(color="black"),
    axis.ticks = element_line(color="black"))


llabs <- c("0" = "Rel. Discrim. = 0",
           "0.333333343267441" = "Rel. Discrim. = .3",
           "0.5" = "Rel. Discrim. = .5",
           "0.666666686534882" = "Rel. Discrim. = .66",
           "1" = "Rel. Discrim. = 1")

mydf <- ggpredict(interactionfullcontrols20, terms = c("rr", "reldep[0, 0.5, 1]"), ci.lvl = 0.95)
int20 <- ggplot(mydf, aes(x, predicted, group = group)) +  geom_line() +
  geom_line(aes(y = conf.high), linetype = 2) +
  geom_line(aes(y = conf.low), linetype = 2) + facet_wrap(~group, nrow = 2, labeller=labeller(group = llabs)) +
  xlab("Racial Resentment") + ylab("pr(Voted)") +
  scale_x_continuous(breaks = c(0,
                                .5,
                                1)) +
  ylim(c(0,1)) +
  ggtitle("Rel. Discrimination *\nRacial Resentment") +
  theme_classic() +
  theme(
    plot.title=element_text(hjust=.5, size=12, face="bold"),
    axis.text = element_text(color="black"),
    axis.ticks = element_line(color="black"))

int20


MFX20 <- (int20 | (rr20 / rd20)) + plot_layout()
MFX20

png("Figure 3.png", width = 750, height = 500, units = "px")
wrap_elements(MFX20) 
dev.off()

ggsave(wrap_elements(MFX20) ,
       filename = "Figure 3.pdf",
       device=cairo_pdf,
       dpi = 300,
       units = "in",
       width = 6,
       height = 4)

##############Rural Consciousness#####################

interactionfullcontrols20 <- glm(voted ~ reldep * rr + age + income + education + female + resimobi + ownhome + strpty + church_attoft + union + contacted , data = whites20, family = "binomial")
summary(interactionfullcontrols20)

interactionfullcontrols20RC <- glm(voted ~ reldep * rr + ruralconsciousness + age + income + education + female + resimobi + ownhome + strpty + church_attoft + union + contacted , data = whites20, family = "binomial")
summary(interactionfullcontrols20RC)

eb_interactionfullcontrols20 <- lm(eb_index7 ~ reldep * rr + age + income + education + female + resimobi + ownhome + strpty + church_attoft + union + contacted, data = whites20)
summary(eb_interactionfullcontrols20)

eb_interactionfullcontrols20RC <- lm(eb_index7 ~ reldep * rr + ruralconsciousness + age + income + education + female + resimobi + ownhome + strpty + church_attoft + union + contacted, data = whites20)
summary(eb_interactionfullcontrols20RC)

ce_interactionfullcontrols20 <- lm(ce_index ~ reldep * rr + age + income + education + female + resimobi + ownhome + strpty + church_attoft + union + contacted, data = whites20)
summary(ce_interactionfullcontrols20)

ce_interactionfullcontrols20RC <- lm(ce_index ~ reldep * rr + ruralconsciousness + age + income + education + female + resimobi + ownhome + strpty + church_attoft + union + contacted, data = whites20)
summary(ce_interactionfullcontrols20RC)

stargazer(interactionfullcontrols20,
          interactionfullcontrols20RC,
          eb_interactionfullcontrols20,
          eb_interactionfullcontrols20RC,
          ce_interactionfullcontrols20,
          ce_interactionfullcontrols20RC,
          title="Table C1: Effects discussed above not conditional on Rural Consciousness",
          align=TRUE,
          dep.var.labels.include=FALSE,
          dep.var.caption = "",
          omit.stat = c("ser"),
          order = c("reldep:rr", "rr", "reldep", "ruralconsciousness", "age", "income", "education", "female", "resimobi", "ownhome", "strpty", "church_attoft", "union", "contacted", "Constant"),
          covariate.labels = c("Interaction Term",
                               "Rel. Discrimination",
                               "Racial Resentment",
                               "Rural Consciousness",
                               "Age",
                               "Income",
                               "Education",
                               "Female",
                               "Resid. Mobil.",
                               "Home Owner",
                               "PID Strength",
                               "Freq. Church Att.",
                               "Union Mem. in HH",
                               "Contacted",
                               "Constant"),
          star.char = c("*"),
          star.cutoffs = c(0.05),
          notes = c("* p<0.05"), 
          notes.append = F,
          out = "Table C1.html")

